• variables: math ability, socio-economic status

• variables: math ability, socio-economic status
• N = 7185 pupils across 160 schools

• what holds true on the group level need not be true on the individual level

• what holds true on the group-level need not be true on the individual level
• need to take both into account

• what holds true on the group-level need not be true on the individual level
• need to take both into account
• leads to hierarchical modeling

\[Y_{i} = \beta_0 + \beta_1 cses_i + \varepsilon_{i}\]

        Mean      SD 
b0    12.762 0.08003 
b1     2.192 0.12087 
sigma  6.725 0.05648 

\[Y_{ij} = \beta_0 + \beta_1 cses_i + \nu_{0i} + \varepsilon_{ij}\]

               Mean        SD 
b0[1]      9.955816 0.8507069 
b0[2]     13.372271 1.0925979 
b0[3]      8.076121 0.8420225 
b0[4]     15.557423 1.2288069 
b0[5]     13.107298 0.8414050 
b0[6]     11.427264 1.0153226 
b0[7]     10.140683 1.0704554 
b0[8]     18.973798 0.9577321 
...
b1         2.192245 0.10753146 
mu.int    12.656230 0.24356510 
sigma      6.085708 0.05064550 
sigma.int  0.117515 0.01433764 

\[Y_{ij} = \beta_0 + \beta_1 cses_i + \nu_{0i} + \nu_{1i} + \varepsilon_{ij}\]

                  Mean        SD 
b0[1]         9.986826 0.8547047 
b0[2]        13.422617 1.1352827 
...
b1[158]      2.4861617 0.75520371 
b1[159]      2.5727856 0.72138232 
b1[160]      2.0203291 0.70894152 
mu.int      12.6357940 0.24790906 
mu.slope     2.1753626 0.12624520 
sigma        6.0597152 0.05286105 
sigma.int    2.9698435 0.19077038 
sigma.slope  0.8257171 0.16246617 

• computationally more feasible

• computationally more feasible
• Bayesian analog to AIC designed for hierarchical models

• computationally more feasible
• Bayesian analog to AIC designed for hierarchical models
• Accounts for flexibility of priors: penalizes more flexible models more than those with more constrained

• computationally more feasible
• Bayesian analog to AIC designed for hierarchical models
• Accounts for flexibility of priors: penalizes more flexible models more than those with more constrained
• tends to favour the more complex model (like AIC; also inconsistent)

Model 1: \(Y_{i} = \beta_0 + \beta_1 cses_i + \varepsilon_{i}\)

Mean deviance:  47777 
penalty 3.134 
Penalized deviance: 47781 

Model 1: \(Y_{i} = \beta_0 + \beta_1 cses_i + \varepsilon_{i}\)

Mean deviance:  47777 
penalty 3.134 
Penalized deviance: 47781 

Model 2: \(Y_{ij} = \beta_0 + \beta_1 cses_i + \nu_{0i} + \varepsilon_{ij}\)

Mean deviance:  46339 
penalty 147.7 
Penalized deviance: 46486 

Model 1: \(Y_{i} = \beta_0 + \beta_1 cses_i + \varepsilon_{i}\)

Mean deviance:  47777 
penalty 3.134 
Penalized deviance: 47781 

Model 2: \(Y_{ij} = \beta_0 + \beta_1 cses_i + \nu_{0i} + \varepsilon_{ij}\)

Mean deviance:  46339 
penalty 147.7 
Penalized deviance: 46486 

Model 3: \(Y_{ij} = \beta_0 + \beta_1 cses_i + \nu_{0i} + \nu_{1i} + \varepsilon_{ij}\)

Mean deviance:  46279 
penalty 190.7 
Penalized deviance: 46470 

• use the same formula for AIC weights (see Wagenmakers & Farrell (2004)), but without the \(\frac{1}{2}\) in the exponent

• use the same formula for AIC weights (see Wagenmakers & Farrell (2004)), but without the \(\frac{1}{2}\) in the exponent
• \(\vartriangle_i (DIC) = DIC_i - \text{min} DIC\)

• use the same formula for AIC weights (see Wagenmakers & Farrell (2004)), but without the \(\frac{1}{2}\) in the exponent
• \(\vartriangle_i (DIC) = DIC_i - \text{min} DIC\)
• \(w_i(DIC) = \frac{e^{-\vartriangle_i (DIC)}}{\sum\limits_{k=1}^K e^{-\vartriangle_{k}(DIC)}}\)

• use the same formula for AIC weights (see Wagenmakers & Farrell (2004)), but without the \(\frac{1}{2}\) in the exponent
• \(\vartriangle_i (DIC) = DIC_i - \text{min} DIC\)
• \(w_i(DIC) = \frac{e^{-\vartriangle_i (DIC)}}{\sum\limits_{k=1}^K e^{-\vartriangle_{k}(DIC)}}\)
• interpretation: model is the best model given the data and the set of candidate models

dic_weights(dic_samples, which = c('rint', 'rslope'))

Inf

dic_weights(dic_samples, which = c('maximal', 'rint'))

33488160

• fit the model using lme4

• fit the model using lme4
• approximate the Bayes factor via BIC (Wagenmakers, 2007)

• fit the model using lme4
• approximate the Bayes factor via BIC (Wagenmakers, 2007)
• \(\text{BIC} = -2 \log p(\theta|y) + k \log n\)

• fit the model using lme4
• approximate the Bayes factor via BIC (Wagenmakers, 2007)
• \(\text{BIC} = -2 \log p(\theta|y) + k \log n\)
• \(BF_{01} = \exp{\frac{BIC_1 - BIC_0}{2}}\)

BF01_bic(freq$rslope, freq$rint)

3.243075e+228

BF01_bic(freq$maximal, freq$rint)

64.27706 # rint over maximal (!!)

BF01_bic(freq$maximal, freq$rint)

64.27706 # rint over maximal (!!)

• apparently the BIC Bayes factor approximation yields the opposite result as the DICw

BF01_bic(freq$rint, freq$maximal)

64.27706 # rint over maximal (!!)

• apparently the BIC Bayes factor approximation yields the opposite result as the DICw
• a classical test is inline with the DICw

Data: dat
Models:
freq$rint: mathach ~ cses + (1 | school)
freq$maximal: mathach ~ cses + (cses | school)
             Df   AIC   BIC logLik deviance  Chisq Chi Df Pr(>Chisq)   
freq$rint     4 46728 46756 -23360    46720                            
freq$maximal  6 46723 46764 -23356    46711 9.4331      2   0.008946 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Data: dat
Models:
freq$rint: mathach ~ cses + (1 | school)
freq$maximal: mathach ~ cses + (cses | school)
             Df   AIC   BIC logLik deviance  Chisq Chi Df Pr(>Chisq)   
freq$rint     4 46728 46756 -23360    46720                            
freq$maximal  6 46723 46764 -23356    46711 9.4331      2   0.008946 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

• sorry for the p-value!

Data: dat
Models:
freq$rint: mathach ~ cses + (1 | school)
freq$maximal: mathach ~ cses + (cses | school)
             Df   AIC   BIC logLik deviance  Chisq Chi Df Pr(>Chisq)   
freq$rint     4 46728 46756 -23360    46720                            
freq$maximal  6 46723 46764 -23356    46711 9.4331      2   0.008946 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

• sorry for the p-value!
• it pains me at least as much as it pains you ;)

\(Y_{ij} = \beta_0 + \beta_1 cses_i + \nu_{0i} + \nu_{1i} + \varepsilon_{ij}\)

                                   95% Highest Density Interval
              Mean        SD         lower     higher
b0[1]         9.986826 0.8547047 
b0[2]        13.422617 1.1352827 
b0[3]         8.035897 0.8398268 
...
b1[158]      2.4861617 0.75520371 
b1[159]      2.5727856 0.72138232 
b1[160]      2.0203291 0.70894152 

mu.int      12.6357940 0.24790906 12.17682  13.11603
mu.slope     2.1753626 0.12624520 1.948548  2.442388
sigma        6.0597152 0.05286105 5.956384  6.149032 
sigma.int    2.9698435 0.19077038 2.626246  3.268999
sigma.slope  0.8257171 0.16246617 0.5430283 1.127743

• CSES is positively related to mathematical achievement
• This relationship differs across schools

• use Savage-Dickey density ratio (more on that later)

• use Savage-Dickey density ratio (more on that later)
• test if catholic schools differ from public schools

• use Savage-Dickey density ratio (more on that later)
• test if catholic schools differ from public schools

BF01_bic(freq$maximal, freq$maxsector)

200507890

savage_dickey(sector_delta, prior = dcauchy(0))

60719.16

• BIC Bayes factor approximation is weird for hierarchical models

• BIC Bayes factor approximation is weird for hierarchical models
• DIC is like AIC (favours complex models, inconsistent)

• trials are nested within participants

• trials are nested within participants
• experimental effect probably isn't the same for every participant

• 2 conditions, within-subject

• 2 conditions, within-subject
• \(\alpha_i\) = \(\phi_{exp}\) - \(\phi_{base}\)

• \(\phi_{base} \sim N(\mu_{\phi}, \sigma_{\phi}^2)\)

• \(\phi_{base} \sim N(\mu_{\phi}, \sigma_{\phi}^2)\)
• \(\alpha_{i} \sim N(\mu_{\alpha}, \sigma_{\alpha}^2)\)

• \(\phi_{base} \sim N(\mu_{\phi}, \sigma_{\phi}^2)\)
• \(\alpha_{i} \sim N(\mu_{\alpha}, \sigma_{\alpha}^2)\)
• \(\mu_{\phi}\), \(\sigma_{\phi}\), \(\sigma_{\alpha}\) are given uninformative priors

• \(\phi_{base} \sim N(\mu_{\phi}, \sigma_{\phi}^2)\)
• \(\alpha_{i} \sim N(\mu_{\alpha}, \sigma_{\alpha}^2)\)
• \(\mu_{\phi}\), \(\sigma_{\phi}\), \(\sigma_{\alpha}\) are given uninformative priors
• \(\mu_{\alpha} = \delta \times \sigma_{\alpha}\)

• \(\delta = \frac{\mu_{\alpha}}{\sigma_{\alpha}}\)
• \(\delta \sim \mathrm{Cauchy}(s = \frac{\sqrt{2}}{2})\)
• see e.g. Rouder & Morey (2012)

library('rjags') # interface with JAGS

library('coda') # MCMC diagnostics

library('ggmcmc') # nice plots

library('polspline') # get height of posterior

• for details, see accompanying R code

ggs_traceplot(ggs(samples))

gelman.diag(samples)

ggs_density(ggs(samples))

• we are not interested in estimating things that may or may not exist

• we are not interested in estimating things that may or may not exist
• parameter testing logically precedes parameter estimation

\(BF_{01} = \frac{p(D|H_0)}{p(D|H_1)} = \frac{p(\delta \, = \, 0|D, H_1)}{p(\delta \, = \, 0|H_1)}\)

library('polspline')

savage_dickey <- function(prior) {
  fit.posterior <- logspline(samples, lbound = 0)
  posterior <- dlogspline(0, fit.posterior)
  BF01 <- posterior / prior
  1 / BF01
}

1 / savage_dickey(prior = dcauchy(0))

3.56

savage_dickey(prior = dunif(1, 0, 1000000))

Inf

savage_dickey(prior = dunif(1, 0, 1000000))

Inf

• Lindley's Paradox (DeGroot, 1982; Lindley, 1957)

savage_dickey(prior = dunif(1, 0, 1000000))

Inf

• Lindley's Paradox (DeGroot, 1982; Lindley, 1957)
• diffuse priors should not be used for hypothesis testing since they lead to unbounded support for \(H_0\)

• simultanously model group- and individual-level effects
• shrinkage
• relaxing assumptions (unbalanced designs, dependencies, etc.)
• many, many more; barely scratching the surface

• natural to think about the group-level in terms of hyperpriors
• just ''turning the handle'' (no REML / ML etc.)

• code is in the ttest/ directory, adapted from Wagenmakers, Lodewyckx, Kuriyal, & Grasman (2010)
• since we work with binary data, we need to model it accordingly
• that is, model the probit of the response
• using ANOVA (with transformations) would lead to spurious results (Jaeger, 2008)

• use titanic survivor data and logistic regression
• compare simulated uniform covariate with Savage-Dickey and BIC Bayes factor
• see comparison/sanity

• compare random effects with DIC weights and BIC Bayes factor
• additionally compare fixed effects with Savage-Dickey
• see comparison/hierarchical/

• in lme4 n is always the number of data points
• however, each level has it's own n
• e.g. 160 schools. test on this level has n = 160
• e.g. 7185 pupils, test on this level has n = 7185
• also, \(k\) is not a good penalty for model complexity